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ABSTRACT 


This contract final technical report presents the theoretical derivation of a finite element solution 
algorithm for the transient energy conservation equation in multidimensional, stationary multi-media 
continua with irregular solution domain closure. It establishes the complete finite element matrix forms 
for arbitrarily irregular discretizations using natural coordinate function representations. It describes 
embodiment of the algorithm into a user-oriented computer program (COMOC) which obtains transient 
temperature distributions at the node points of the finite element discretization using a highly stable ex- 
plicit integration procedure with automatic error control features. The finite element algorithm is shown 
to possess convergence with discretization for a transient sample problem. The condensed form for the 
specific heat element matrix is shown to be preferable to the consistent form. Sample results computed 
for a practical aerospace problem indicate excellent solution accuracy and speed of computation, and 
factors within the integration package which affect its operation are evaluated. Computed results for 
diverse problems illustrate the versatility of COMOC, and easily prepared output subroutines are shown 
to allow quick engineering assessment of solution behavior. 
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I. INTRODUCTION 


Numerical solution of the discretized equivalent of the initial bouhdary value problem of mathemati- 
cal physics is commonplace with the advent of the large digital computer. The finite element method, well 
proven for solution of large-scale problems in structural mechanics, is being extended to general boundary 
value problems throughout continuum mechanics. For these problems in heat transfer, fluid mechanics, 
electromagnetic theory, and magneto hydrodynamics, application of the finite element solution algorithm to 
the original system of quasilinear partial differential equations produces large-order systems of ordinary dif- 
ferential equations written on the discretized equivalent of the original dependent variables. These equations 
are eligible for solution by both implicit and explicit numerical integration techniques. 

Recent research has yielded analytical techniques for deriving integration methods for ordinary dif- 
ferential equations that are optimum on the multiple bases of stability, accuracy and required computing time. 
Computational methods are available for establishing specifically tailored 1-, 2-, and 3-step explicit integra- 
tion algorithms, satisfying both accuracy and stability requirements if possible. These techniques have been 
used to establish a class of explicit numerical integration methods which are one-step (and therefore self- 
starting), have error control features, require derivative evaluations at the integration-interval end points only, 
are optimally stable and accurate within their given structure, and are readily programmed for use on the 
digital computer. These methods represent a generalization of the well known one-step predictor-corrector 
methods with 2 corrections and 3 derivative evaluations per step. 

The results of numerical evaluation of the highly stable integration algorithm are presented herein as 
established via solution of the oridinary differential equation system resulting from a simplex finite element 
discretization of two-dimensional and axisymmetric problems in transient heat transfer. These results cor- 
roborate the theoretical preditions and illustrate the usefulness of the automatic error control features. The 
versatility of the finite element solution technique is similarly exemplified for practrical problems involving 
complex geometries, tensor conductivities, and temperature dependent thermophysical properties. 
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II. FINITE ELEMENT SOLUTION ALGORITHM FOR 
ENERGY CONSERVATION 


A. GOVERNING DIFFERENTIAL EQUATIONS 

The fundamental concept in continuum mechanics analysis is identification of material distributions 
over domains, and satisfication of the conservation laws of mechanics expressed in terms of the continuum 
distributions. For a bounded domain at rest in an inertial reference frame, containing material of constant 
but non-uniform density, the conservation laws not identically satisfied reduce to the first and second laws of 
thermodynamics. In its elementary form, the first law of thermodynamics states 

AE = AQ- AW (2.1) 

i 

where AE is the energy change within a domain when AQ amounted heat has been added and AW amount of 
work has been done. Averaging Eq (2.1) over a time interval At and taking the limit as At approaches zero, 
we have, assuming the limit exists and is finite, the equivalent differential equation statement of the first law. 

dE cTQ dW 

= AAA - I (2.2) 

dt dt dt 


In the context of the present development, the second law of thermodynamics places a lower limit on allowed 
entropy changes between thermodynamic end points. 


Eq (2.2) is valid over all domains. It can be specialized to a continuum distribution by identification 
of appropriate thermodynamic and material variables. The internal energy E is expressed in terms of a cor- 
responding internal energy density, e, as 



The thermodynamic point function, e, can be expressed as a function of two other variables, say density and 
temperature. Since density is at least piecewise constant, and identifying a specific heat capacity, c, obtain 
for the internal energy density. 


e 


J c(T) dT 


(2.4) 


Considering the remaining terms in Eq (2.2), the work done by an incompressible continuum at rest 
is zero. The heat-added term may have contributions stemming from internal generation, Q, and a net inward 
flux of heat, q^, over the bounding surface dR of the domain R. Concerning the latter, it is convenient to 
relate forces and resultant fluxes by material coefficients. Accounting for the vector character of heat flux 
and temperature gradient, the material property flux coefficient can at most be a second-order tensor. Using 
Cartesian tensor subscript notation, with the summation convention implied for repeated Latin indices and 
partial differentiation denoted by a comma, the generalized flux form is 

<U = k ij T >j (2.5) 

where k^, the thermal conductivity tensor, is typically temperature -dependent for practical problems. 
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Employing the defined material properties, and the concept of mass distribution, the first law of 
thermodynamics, Eq (2.2), can be written as an integro-differential equation of the form 



nj da 




p cdTdr 


(2.6) 


Employing Gauss’ theorem on the surface integral, and noting that the domain R is time independent and 
arbitrary, the integrand must vanish identically. The familiar differential equation for energy conservation 
at a point is then 

(ky T ;j ), j + pQ = pcT, t (2.7) 

With the explicit assumption that material properties can be discontinuous, Eq (2.7) is valid throughout R 
by subdomains with boundaries identical with the surfaces of these discontinuities. 

Equation (2.7) is an elliptic boundary value problem that is also initial value. Hence, a unique solu- 
tion exists for appropriate specification of boundary and initial conditions. The initial condition is simply 

T (Xj, 0) = T q (Xj) (2.8) 

The generalized elliptic boundary condition is a first order differential equation relating the temperature 
and its normal derivative, on the complete closure of the domain R, of the form 

aj (t) T (x i; t) + a 2 (t)T(xj, t), k n k = a 3 (t) (2.9) 

In Eq (2.9), should a, vanish on any closure segment, a fixed heat flux is applied to that portion of the do- 
main closure. The vanishing of a 2 enforces a fixed temperature. Should none of the aj vanish, the temperature 
and its normal gradient are linearly related, the so-called convection boundary condition, usually written as 

% (xj, t) = h (Xj, t) [ T (Xj, t) - T r (xj, t)J nj (2.10) 

where h is the defined convective heat transfer coefficient and T r is the temperature of the reference exchange 
medium. 


Another form of surface energy flux with the domain of interest, requiring no close communication 
with the source, is radiant energy exchange. This boundary condition statement is similar to Eq (2.10), but 
with fourth power appearance of the temperature terms and generalization of the heat transfer coefficient to 
include geometric viewfactors. The temperature quartic can be factored once to yield Eq (2.10) explicitly, 
with the following expression resulting for the radiation heat transfer coefficient. 

h r (Xj,t) = a a (Xj, t) [ T 3 + T 2 T f + TT r 2 + T f 3 ] (2.11) 

In Eq (2.1 1), a is the Stefan - Boltzmann constant and a is the viewfactor. Eq (2.9) is thus noted to embody 
both of the familiar convection and radiation boundary condition forms. 

B. FINITE ELEMENT SOLUTION ALGORITHM 

For a numerical solution algorithm of the equation system (2.7) - (2.9) to be practical, it must 
explicitly admit variable thermophysical properties and arbitrarily shaped solution domain closures. The 
finite element procedure is distinctly appropriate for this purpose. Use is made of Galerkin criteria in the 
Method of Weighted Residuals (MWR), on a local basis, to provide the theoretical foundation for the finite 
element solution algorithm. 
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Equation (2.7), a quasilinear elliptic partial differential equation, can be written in the form 


L(T) = (kjj T, j), j + f(T,T,j 


Xj,t) = 0 


( 2 . 12 ) 


Classical application of Galerkin criteria within MWR constitutes forming an approximation, T* (Xj, t) to the 
dependent variable, T (Xj, t), by a series expansion into a set of unknown coefficients modifying members, 
0 k > of a complete set of functions in the spatial variables. i l 


K 

T* (*i> t) = S C k (t) 0 k (x A ) 


(2.13) 


Certain of these expansion coefficients are evaluated by requiring that the approximate solution satisfy the 
boundary conditions, Eq. (2.9). The remaining unknown coefficients are then determined by setting to zero 
the differential equation residuals, formed by substitution of T* into Eq (2.12), multiplication by a set of 
weighting functions, W k , typically identical to the spatial approximation functions, 0 k , and integrating over 
the domain R. This produces the N ordinary differential (or algebraic, for a time independent problem) 
equations . I 


i 


W k (x i )L(T*)dr=0 k=l[.2, ...,N 


(2.14) 


where N corresponds to the total number of remaining unknown coefficients in the definition of the approxi- 
mation function, Eq (2.13). 


When non-vanishing gradient boundary conditions exist, a generalization to the classical MWR ap- 
proach is advantageous, to relax the constraints formed by the boundary condition statement, Eq (2.9). 
Form the N boundary residuals, after generalization of a 2 , Eq. (2.9), as 


Lr 


w k ( x i) [T* (a 3 - a 2i j T, j nj)] da = 0 k = 1 ,2, , . N 


(2.15) 


Looking to the variational calculus for guidance, multiply Eq (2.15) by a Lagrange multiplier, X, and sub- 
tract them from Eq. (2.14). Noting that the N weighting functions can be written as, 


W k (XJ) = 


3T* ( Xj , t) 

ac k (t) 


(2.16) 


integrating by parts the differential equation term involving the elliptic operator, identifying X with unity, 
and specifying that a 2 coincide with the thermal conductivity tensor, a cancellation of terms is achieved. The 
resultant equation system for determination of the remaining expansion coefficients in Eq. (2.13) becomes 



f (t* j t*, j, Xj, t) dr - 



9T* 

Ui T* -a 3 ] do = 0 

dC k 


k = 1, 2, . . . , N (2.17) 


It is readily shown that Eq (2.17) is identical to that equation formed by extremization of the variational 
equivalent statement of Eq (2.7) - (2.9), when the equations are linear and f is explicitly independent of 
both T and t. 
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The Numerical Method of Weighted Residuals follows directly from this theory. Form a series ap- 
proximation, T* m , to the dependent variable, T, constrained to lie within the m th subdomain, R m , defined 
by (x-, t) e R x f 0, °°) , and where the union of the M disjoint interior subdomains, Rp, forms R. Evaluate 
Eq (2T7) in each of these “finite elements”, R m , applying the gradient boundary condition statement to 
those elements with surfaces coincident with the solution domain closure 3R. Assemble the M x N Eq 
(2.17) into a global system written on the totality of unconstrained expansion coefficients, Cp k (t), by 
Boolean algebra. The equation system thus produced is again identical to that formed by a variational pro- 
cedure when a stationary principle exists for Eq (2.7) - (2.9), i.e., when the differential equation system is 
linear and time-independent. HQwever, since no linearity constraint has been employed in the derivation of 
the solution theory, it is generally applicable to all problem descriptions belonging to Eq (2.7) - (2.9). 

The final step for establishment of the finite element solution algorithm is selection of the func- 
tionals, <£ k (Xj), Eq (2.13), and determination of the desired shape of the discretized subdomains, R m - The vast 
experience in structural analysis has yielded finite element shapes and approximation functions which are 
highly specialized to each particular problem class. At this juncture, for the energy conservation problem, it 
appears that use of low order polynomials spanning the space of the subdomain can provide sufficiently ac- 
curate and economical solutions for practical engineering problems. 


The further restriction that 0 k (xj) be linear polynomials establishes the finite element shapes shown 
in Figure 1 for 1-, 2- and 3-dimensional spaces respectively. Selecting node locations on the closure of the 
subdomain, 3R m , as indicated, allows physical significance to be attached to the expansion coefficients, Eq 
(2.13), by requiring that the series represent (dimensionally) the temperature at each node. With this special- 
ization, the approximation to temperature in the m^ finite element takes the form 

IT 


T*m ( x i> t) = L 


e (t) 


m 


(2.18) 


The elements of the vector j 0 (t)[ are the time dependent temperatures at nodes of the discretization and 
the elements of j L| are linear polynomials in the spatial coordinates. With identification of Eq (2.1 8) the 
weighting functions, Eq (2. 1 6), for Galerkin criteria are simply the elements of | L } . Combining Eq ( 1 .27) - 
(2.18), the finite element solution algorithm for Eq (2.7) - (2.9), within the finite element R m , can be 
compactly written as 


[C] 


m 


{#(t)} 


m 


= - 2 

I 


lKIl m {<><t)} m + S {Yl} 


m 


(2.19) 


by identification of the following element matrices. 
[CJ m h 


[Kl] m 

[K2) m 

{Yl} m 

|Y2} 


/ l L f Pm c m l L f T dT 

' R m 

[ |l} , jfcqj 171 ] |l| T ,j dr 

J R m 

/ 


JR m 


l L l a i m W da 

W p m Q m dr 


m 


I l L f a 3m da 


( 2 . 20 ) 
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In Eq (2.19) the superscript prime denotes the ordinary time derivative of the element vector of nodal tem- 
peratures. The finite element matrices defined in Eq (2.20) are respectively the specific heat capacity, the 
thermal conductivity, the temperature dependent contribution to the flux boundary condition, the internal 
heat generation, and the temperature independent contribution to the flux boundary condition. 

Equations of the form (2.19) are written for every finite element lying within the solution domain. 
Assembly by Boolean algebra of the M Eq (2.19), into the required global solution form, yields an equation 
system which is identical in appearance, save for deletion of the subscript m. Global equation (2.19) is 
formally rendered in standard form, for solution by any (finite difference) integration algorithm, by pre- 
multiplication with [Cj -1 , the inverse of the global heat capacity matrix. Thus the final form of the finite 
element solution algorithm for the time dependent distribution of node temperatures of the discretization 
becomes 

|0(t)f' = [Cj" 1 (- 2 [KI] {0(t)| + 2 |YI|) (2-21) 

Boundary condition specifications are accepted on an element basis by the set of aj m , Eq (2.9), for 
each finite element with nodes occurring on the closure of the global solution domain. For a typical ele- 
ment, should a t m be zero, the boundary temperature gradient is prescribed (fixed heat flux) and accounted 
for by { Y2} , Eq (2.20). Should a 2 m vanish identically, the node temperature is fixed and not calculable 
from Eq (2.21). Finally, if none of the aj m vanish, the contributions are contained within both the { Y2 } 
and [K2] matrices, Eq (2.20). The more familiar convective boundary condition statement, Eq (2.10), 
results for a, m identified with the heat transfer coefficient, h, and a 3rn set equal to a, m multiplied by the 
reference temperature for heat exchange. Thus, the finite element algorithm accepts an arbitrarily varying, 
well-posed elliptic boundary condition statement on segments of all finite elements coincident with the 
global solution domain closure. 

C. FINITE ELEMENT MATRIX GENERATION 

Implementation of the solution algorithm into an operational computer code basically involves evalu- 
ation of the finite element matrices, Eq (2.20) and assembly and integration of Eq (2.21 ). For the simplex 
1-, 2-, and 3-dimensional elements chosen, Figure 1, the use of natural coordinate systems is preferred, since 
evaluation of the required moment distributions, Eq (2.20), is particularly straightforward. A simplex natural 
coordinate is a (linear) function spanning the domain of a finite element, with the requirement that it vanish 
at all nodes of the domain except one where its value is unity. These functions are discussed in some detail 
in Reference 1, and are generalizations of the area coordinates of structural mechanics, extended to 1-, 2- and 
3- dimensions, see Zienkiewicz 

The definitions of the natural coordinates for simplex functionals spanning 1-, 2- and 3- dimensional 
spaces are expressed in Table 1 . The { L} coordinate system is defined implicitly in terms of linear series 
approximations of the function involving the node point coordinates of the finite element. In each case, the 
first definition equation enforces the linear dependence of the natural system. The determinant of the 
matrix of coefficients modifying the definition of is algebraically related to the length, plane area, and 
volume of the 1-, 2- and 3-dimensional elements respectively. Evaluation of the finite element matrices, 

Eq (2.20) uniformly involves moment distributions over the domain of the finite element. The values taken 
for these moments are listed in Table 2 and are noted to involve only the exponent distribution of the ele- 
ments of { L} . Note also the commutativity of the exponent distribution. 

1 . Planar Elements in Two-Dimensional Space 

Evaluation of the finite element matrices in 2-dimensional space is straightforward, since mo- 
ment generation is invariant under coordinate translation and rotation. The preferred orientation for the 
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TABLE 1 

IMPLICIT DEFINITION OF LINEAR NATURAL 
COORDINATE FUNCTIONS 



TABLE 2 

INTEGRALS OF THE NATURAL COORDINATE FUNCTIONS 


Dimensions 

Integrals* 


f 01 

n 2 


1 

/ Ll 

l 2 

dx 


J* 




f 

n 2 

n 3 

2 

/ L ' 

1-2 

L 3 dxdy 


Jr 




f n, 

n 2 

n 3 n 4 

3 

1 L ' 

La 

L 3 L 4 dxdydz 


R 




= D 


nj ! n 2 ! 
(1+rij + n 2 )j 


D 


n i I n 2 ! n 3 ! 

(2+ nj +n 2 + n 3 )! 


n i ! n 2 ! n 3 ! n 4 ! 

(3+ n x +n 2 +n 3 +n 4 )! 


* D = Determinant of coefficient matrix defining the 
natural coordinate system, see Table 1. 
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triangular element is in a local (unprimed) coordinate system, Figure 1 , defined with respect to the global 
system by the tensor transformation law 


Xj = ^ x ; 

The volume integration kernel, Eq (2.20), becomes 


( 2 . 22 ) 


dr = t dxdy (2.23) 

where t is the thickness of the finite element. Surface integrals may be of two types involving integration 
over the plane area or along an edge of an element. In these cases, the surface integration kernels, do, 
respectively take the forms 

! dxdy 

(2.24) 

tdx 

where the x axis is always chosen to coincide with the side of the element coincident with the global domain 
closure segment. Except for the thermal conductivity matrix, [ K 1 ] , the evaluations are routine. Shown 
in the first column of Table 3 are the finite element matrices for all of Eq (2.20), including the element 
dependent algebraic multiplier. A m is the plane area of the finite element, and the boundary condition 
matrices are specified in terms of the convection coefficients h m and T^ and the specified heat flux F m . 


The thermal conductivity matrix requires some additional comment. Rewriting the second 
of Eq. 2.20 for 2-dimensional space, the [K1 ] matrix becomes 


[Kl] m 



i L }>j kjj™ {L} (i t dxdy 


(2.25) 


where ky is the element thermal conductivity tensor defined in the unprimed coordinate system. Since 
the elements of |l[ are linear, their first derivatives are constants and the integration reduces to evalua- 
tion of the plane area of the finite element. Performing the indicated summation operations on the repeated 
Latin subscripts, the [K1 ] matrix takes the final form. 



L l,l 

r 

<N 

s 


m 
ki i 

i 

E « 

[Kl] m J ^m t m 

L 2,l 

L 2,2 


m 
k 2 1 

m 
k 2 2 


- L 3,l 

L 3,2_ 





L 1 ,1 L 1 ,2 L 1 ,3 

_ L 2,1 l 2,2 l 2,3_ 


(2.26) 


The ky tensor scaler components are obtained from the transformation definition and Eq (2.22) as 


m m' 

ky - api a^j kpq 


(2.27) 


Therefore, evaluation of [K1 ] involves differentiation of the natural coordinate functions, Table 2, and 
the indicated matrix multiplications. 
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TABLE 3 

FINITE ELEMENT MATRICES FOR THE ENERGY CONSERVATION EQUATION 
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2. Planar Ring Elements in Axisymmetric Two-Dimensional Space 

In this space, the matrix forms, Eq (2.20), are invariant only under coordinate rotation, and 
require redefinition of the integration kernels. From the Jacobian of the transformation, Eq (2.22), and 
differential geometry, the integration kernels in the unprimed coordinate system become 

dr = 2irr' dxdy (2.28) 

da = 27rr' dl 

where r' refers to the radial coordinate of points lying within the finite element expressed in the global 
coordinate system. From the definition of the natural coordinate system, r' becomes 

x =x\ + x, L! + x 2 L 2 +x 3 L 3 (2.29) 


Hence, for example, the specific heat matrix, Equation (2.20) takes the form 

(r{ +x 1 Li+x 2 L 2 +x 3 L 3 ) 27r dxdy (2.30) 


[C] m = P m c 


m m 


Si 

m 


r 2 

L 1 Li L 2 Li L 3 


2 

L 2 


(sym) 


L 2 l 3 

2 

L 3 


Performing the indicated integrations, Table 2, and collecting terms, the final form for [C] m becomes, in 
terms of the nodal radial coordinates (in the global system, but with primes deleted) 



As shown in Figure 1 , node 1 is always the origin of the unprimed coordinate system and the numbering 
is counterclockwise. Note that the first matrix in Eq (2.31) is identical, to within 2 tt, to its 2-dimensional 
counterpart. 

Evaluation of the remaining finite element matrices for axisymmetry, Eq (2.20), is similarly 
straightforward. The final forms are shown in the second column of Table 3. The thermal conductivity 
matrix, which is identical in 2-dimensional and axisymmetric space, is the sole matrix requiring numerical 
matrix multiplications. 
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III. HIGHLY STABLE EXPLICIT NUMERICAL INTEGRATION ALGORITHM 


The explicit numerical solution of stable systems of differential equations with large Lipschitz 
constants creates serious integration step-size restrictions. This has motivated the derivation ot explicit 
methods with increased ranges of stability, see for example References 3-6, as well as a proliferation ot 
implicit methods. Restricting attention to explicit methods, techniques have become available for deriving 
integration methods that are optimum on the multiple bases of stability, accuracy, and required comput- 
ing time( 7 X Methods are now available for evaluating any specified explicit integration algorithm for 
application to solution of a particular system of ordinary differential equations via use of a computer pro- 
gram specifically designed for this purposed). This work was recently specialized to optimally-stable, 
k-stage, one-step numerical integration methods(9), as well as a further restriction to 3-stage, one-step 
methods(lO). The latter methods are distinctly applicable to the ordinary differential equation systems 
resulting from either finite element or finite difference discretizations of elliptic boundary value problems 
that are also initial value. 

Briefly, a generalized family of numerical integration methods that are one-step, predictor- 
multiple-corrector formulas, are described by the equations 

pJi+i = a iy n + hb !yn 

2 2 2 1 2 

P n+1 = a, y n + h[b, p n+ 1 '+-b a y„] 


k - 1 

Pn+1 


k-1 k-1 k-2 , 

a, y n + h[bj p n+1 + 



y n +i 


k-1 


ai y n + h[b,p n+1 '+b 2 



(3.1) 


A member of this family is called a “k-stage, one-step explicit integration method,” and the special case 
k=3 is briefly summarized in the following discussion. The interested reader is referred to Reference 10 
for details. 

The computational structure displayed by Eq (3.1 ) is preferable to the popular k-stage, Runge- 
Kutta scheme, since it contains the potential to approximate the required solution at each stage, thereby 
rendering internal error control a natural adjunct to the process. As with Runge-Kutta methods, 
algorithms of the type Eq (3.1 ) are self-starting and easily programmed for digital computation. The 3- 
stage methods derived in Reference 10 enjoy optimally large stability ranges and minimal truncation error 
bounds. The most stable 3-stage methods are shown to be first-order, and their associated truncation 
errors may be computed with no requirement for additional derivative evaluations. Specialization of Eq 
(3.1) to 3-stage is straightforward, with y interpreted as the numerical approximation to the dependent 
variable, Y, which is the solution to Eq (3.2), 
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= f(Y, t) 


(3.2) 


r dy 
Y = — 
dt 


k k 

In Eq (3.1) li = t n+ j - t n , is the integration step size, and the and tq , k = 1 ,2,3, are real constants. 

The superscript prime denotes evaluation of the numerical approximation of the derivative expression, 
Eq (3.2), and f( ) is a specified function of its argument. Note that Eq (2.2 1) is a vector equation with 
components in the form of Eq (3.2). 


Two members belonging to the 3-stage family, Eq (3.1), have been derived for application to solu- 
tion of large systems of ordinary differential equations. Both methods are first order accurate, i.e., their 
associated truncation error is of order h 2 , and they represent the most stable methods possible within the 
collection of first order accurate methods. The coefficients in Eq (3.1 ) for these two methods are listed 
in Table 4. Method 1 enjoys a large absolute stability interval, while Method 2 has an extended relative 
stability interval. That is, the application of Method 1 for solution of Eq (3.2) provides a propagated 
error which damps to zero for larger values of h than is possible by other methods. Similarly Method 2 
creates a propagated error which damps monotonically for larger values of step size h. The real extent 
of the absolute and relative stability intervals of both integration methods is also shown in Table 4. 


A feature enjoyed by the integration package built around Eq (3.1) is automatic error control 
based upon a user-specified parameter. The relative truncation error estimate magnitude for the subject 
integration methods can be written as ( 1 


I RTE I = - 
(3 


lp n+T -y^l 
IVn+ll 


(3.3) 


The parameter, /3, equals 3 and 6, respectively, for Methods 1 and 2. Eq (3.3) is utilized within the 
integration package to evaluate the relative truncation error associated with using the step size, h, to 
estimate the n+ist value of the dependent variable. If the computed error is less than the user-supplied 
acceptable limit, the n+ist estimate is accepted. At the same time, the integration algorithm will auto- 
matically increase h, by some fraction, before proceeding to the next time-step computation. In this 
fashion, the integration package consistently seeks to increase step-size and hence decrease solution 
computation time. If at some point the computed relative error exceeds the specified limit, the current 
predicted-corrected values for the dependent variable are discarded, a smaller step-size selected, and the 
operations of Eq (3.1) are repeated until an acceptable error is measured. These automatic features may 
be selectively overridden to march at a user-fixed step-size, or at a value of h not exceeding a specified value. 
The user may additionally select an initial step-size to start, or allow the program to compute an acceptable 
initial value. 
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TABLE 4 

COEFFICIENTS IN FINITE DIFFERENCE INTEGRATION ALGORITHM, 
EQ. 3.1 , FOR TWO, ONE-STEP, THREE-STAGE METHODS 


Coefficient 

Method 1 

Method 2 

a i 1 

1.0000 

1.0000 

ai 2 

1.0000 

1.0000 

ai 3 

1.0000 

1.0000 

bi 1 j 

0.1549 

0.2877 

b? 

0.2711 

0.2920 

b 2 2 

0.7289 

0.7080 

bi 3 

0.1667 

0.3334 

b 2 3 

0.8333 

0.6666 


Absolute Stability (-15.997,0] (-8.765,0] 

Interval 

Relative Stability (-1.245,0] (-7.998,0] 

Interval 
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IV. FINITE ELEMENT COMPUTER PROGRAM 


The general purpose finite element computer program COMOC (Computational Continuum 
Mechanics) is coded completely in terms of generalized, nondimensional independent and dependent; 
variables. In addition to producing the results using the Thermal Analysis Variant of COMOC, it has been 
exercised for problems in hydrodynamics (1 O, three-dimensional boundary layer flow02) ; and the 
compressible Navier-Stokes equations^ 3). The program is briefly described as follows; the interested 
reader may obtain details in Reference 14. 

The COMOC program consists of four basic Modules, as illustrated in Figure 2. In INPUT, the 
desired discretization is formed by user specification of the plane global coordinates of the nodes of 
each triangular-shaped finite element. The node numbering sequence is completely arbitrary, as is the 
selection of a particular discretization, and both may be chosen for output convenience. In particular, 
nodes may first be placed where output information is desired. The remainder of the solution domain is 
then discretized to avoid finite elements with aspect ratios greater than about 70:1. The closure of the 
solution domain is recognized by the program as the counterclockwise connection of a specified sequence 
of N nodes. 

The element definition is obtained by identifying the three vertex nodes, the thermophysical 
properties material number, and thickness of each finite element. The non-vanishing boundary condition 
coefficients, aj m , Eq (2.20) are also included on the element cards as is internal heat generation. The 
input phase is completed by reading the material properties tables, various control parameters of the 
integration algorithm of Module 3, initial conditions information, and the desired output control para- 
meters. 


The next two Modules, Figure 2, are basically DO loops on the finite elements of the discretization. 
In the GEOMETRY Module, the finite element matrices, Eq (2.20), are formed. The INTEGRATION 
Module embodies the numerical solution algorithm for both algebraic and ordinary differential equation 
systems; the Thermal Analysis Variant uses only the latter. The basic element operation is formation of 
Eq (2.19), at a given time station, and assembly of the element matrices into the global representation, 

Eq (2.21 ). At user-selected points, the OUTPUT Module is called to record the present values of the 
elements of the temperature vector 1 9 (t)} and other desired parameters. 


COMOC 


Computational Continuum Mechanics 



Figure 2. COMOC Computer Program Organization 
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V. NUMERICAL RESULTS 


An extensive literature has developed over recent years devoted to experimentation with finite 
element methods for problems in heat conduction and other linear field problems. Zienkiewicz and 
Cheung U 5 ) published an early contribution detailing a finite element solution algorithm for a linear 
Poisson equation, and illustrated its usefulness for the steady-state heat conduction problem. Shortly 
thereafter, Wilson and Nickell (16), employing a pseudo-variational principle after Biot ( 17 ) and Gurtin 
established a finite element algorithm for transient heat conduction in linear continua. Baker (* 9 ) 
generalized the theory to include temperature-dependent thermophysical properties. Concurrently, 
Finlayson and Scriven (20, 21) concluded that pseudo-variational principles were generally not required, and 
that comparable algorithm forms resulted from select use of specific criteria within MWR. The 1969 
ASME Winter Annual Meeting Heat Transfer Session (22-27) included many finite element contributions. 
Mason, et.al. (28, 29) employed structural computer codes for thermoelastic analysis of laser mirror 
systems. The use of isoparametric finite elements has been suggested by Zienkiewicz and Parekh ( 3 °) for 
transient heat conduction problems. 

It is quite logical that these referenced authors, whose predominant residence is in the structural 
mechanics community, would contribute to early developments, since they possessed operational finite 
element computer programs that could be adapted for experimentation on heat conduction problems by 
redefinition of dependent variables and continuum properties. These encouraging initial numerical results 
provided the impetus to consider building a computer system tailored to the generalized field problem 
class, and several are now in progress, including the NASTRAN (NASA Structural Analysis) system. This 
is the concept behind COMOC as well. The Thermal Analysis Variant of this computer program has been 
extensively exercised to validate the mathematical soundness of the generalized theory, to evaluate the 
accuracy and convergence character of computed solutions, to prove the significant speed with which 
practical problem solutions can be obtained, and to illustrate the great versatility of the generalized finite 
element computer program system. COMOC has produced the favorable assessments that are discussed 
in this section. 

A. ACCURACY AND CONVERGENCE WITH DISCRETIZATION 

Since no rigorous variational statement exists for the general quasilinear elliptic-initial value problem 
for which the finite element solution algorithm was derived, a numerical evaluation of accuracy of the 
computed solutions, and convergence toward the correct solution with finer discretization, is required. 

The problem selected was external transient convective heat exchange between a medium at constant 
temperature and a thick axisymmetric cylinder with temperature independent thermophysical properties. 

A truncated series solution ( 3 1 ), valid for small solution time and for locations near the exterior cylinder 
surface, was used as the comparison analytic solution. 

The accuracy of solution was found to be a function of discretization as well as the manner select- 
ed for assembly of the global specific heat capacity matrix, [C] . It is recalled that the inverse of [C] is 
required to render the finite element solution algorithm in standard form, Eq (2.21), for integration by 
explicit techniques. Employing the rigorous definition for specific heat capacity matrix, [C] m , Eq (2.20) 
at the element level yields a global [C] which is a sparse symmetric matrix banded about the main diagonal. 
The bandwidth is a function of discretization, and the matrix may be inverted by classical algebraic means, 
or essentially evaluated using equation solver techniques, e.g., banded Cholesky. The [C] matrix may be 
rendered diagonal, and the inversion operation trivial, by condensing the specific heat capacity of each 
finite element at the node points. Such an operation, commonplace in some structural dynamics problems 
wherein lumped mass matrices replace the consistent formulation, is accomplished by adding rowwise 


16 



the matrix elements of each [C] m within a finite element, and placing the nonvanishing coefficient on the 
diagonal. From Table 3, this operation is noted to equally distribute the element heat capacity at the nodes 
for two-dimensional problems; in axisymmetry, the result of condensing is dependent upon the particular 
element location and geometry. 

Both of these options have been evaluated by COMOC for several discretizations of the thick 
cylinder problem. An illustration of the 32 finite element discretization of the cylinder cross-section is 
shown in Figure 3, as well as the various vanishing gradient and convection boundary conditions. Discreti- 
zations employing elements spanning double and half the radial domain (16 and 64 finite elements) of the 
32 element case were also evaluated. The thermophysical material properties were that of Type 304 stain- 
less steel, and the convective heat transfer coefficient and exchange temperature were 34.6 w/m 2 
(20 Btu/hr ft 2 ) and 3 1 1°K (100°F) respectively. The initial cylinder temperature was uniformly 255°K 
(0°F) for all tests. 

Shown in Table 5 is the comparison between the exact solution and the results of computations by 
COMOC, for the transient temperature distribution at and near the surface of the cylinder and within the es- 
sential accuracy bounds of the truncated series solution. The results in the first three COMOC columns, for 
16, 32 and 64 element discretizations using the condensed [C] form, show smooth convergence, generally 
from below, to the correct solution. The same general behavior is noted in the last columns of Table 5 for 
the 16 and 32 element cases using the full [C] form; however, explicit violations of the second law of 
thermodynamics are noted, in both cases, for nodes interior to the surface. At these locations, depression 
of temperature below the uniform 25 5° K (0°F) initial value for a positive heat flux into the solution 
domain is observed. Similar results have occurred for other discretizations and different boundary con- 
ditions and thermophysical properties. In general, a larger applied flux and/or smaller thermal conductivity 
render the temperature depressions more severe (up to -5.0°K (-10°F) for the cylinder made from an in- 
sulator like asbestos). A finer discretization has always improved this adverse phenomena. 

The accuracy comparison between the full and condensed specific heat matrix forms is shown in 
Figure 4. For full [C] , the finite element solution agrees well with the analytic solution at the exterior 
wall, r = 0.61 meters (2.0 ft). The extent of the interior temperature depressions is evident. The 32- 
element, condensed (C ] solution enjoys accuracy comparable to the 16-element full [C] solution at the 
wall, and is uniformly more accurate in the interior. The 16-element condensed [C] solution is considerably 
poorer than either discussed solution. 

These results are in essential agreement with the experience using consistent and lumped mass 
matrix solutions in structural dynamics, as well as the cursory results of Reference 16. The consistent 
(or full) formulation can provide more accurate solutions using coarse discretizations. Both full and con- 
densed solutions converge with discretization; for this particular sample case, a doubling of the discreti- 
zation and moving to the condensed forn>provided a uniformly comparable or more accurate solution. 

These discussed accuracy and convergence results were obtained for the highly stable integration 
algorithm package (QKNINT, QuicK Numerical INTegrator) operating well with its stability bounds, since 
the allowed step-size was constrained by requested output. To determine the stability character, and 
resultant computer CPU execution times for the different [C] matrix options, as well as illustrating the 
automatic step-size and error control features of QKNINT, the convective loading to the 32 element 
thick cylinder was increased 20-fold and the thermal conductivity increased by a factor of 4 to produce 
the required test problem. Figure 5 shows comparison surface temperature computations, and the auto- 
matically determined integration intervals used by QKNINT, Methods 1 and 2 in obtaining the plotted 
transient solutions. Using a value of e = 0.001 for the QKNINT error control parameter, the solution 
accuracy for Methods 1 and 2 for the full [C] matrix form is essentially identical. However, as shown, 
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Figure 3. Discretization of Thick Axisymmetric Cylinder into Thirty-two Triangular Cross-Section 

Ring Finite Elements 
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TABLE 5 

TRANSIENT TEMPERATURE DISTRIBUTIONS IN 
THICK AXISYMMETR1C CYLINDER 


Node 

Radius 

(m) 

Solution 

Time 

(sec) 

Temperature°K 

Series 

Solution 131) 

COMOC Computed Temperature (°K) 

Condensed C Matrix 

Full C Matrix 

16 Element 

32 Element 

64 Element 

16 Element 

32 Element 

0.61 

0.0 

255.38 

255.38 

255.38 

255.38 

255.38 

255.38 


79.2 

259.92 

258.49 

259.58 

259.93 

259.51 

259.85 


158.0 

261.63 

260.65 

261 .50 

261.68 

261.53 

261.66 


237.6 

262.87 

262.15 

262.82 

262.97 

262.84 

262.94 


316.8 

263.88 

263.30 

263.87 

264.02 

263.84 

263.98 


396.0 

264.73 

264.20 

264.75 

264.90 

264.57 

264.87 


475.2 

265.47 

265.20 

265.53 

265.69 

265.31 

265.66 

0.58 

0.0 

255.38 

255.38 

255.38 

255.38 

255.38 

255.38 


79.2 

255.54 

255.53 

255.52 

255.49 

255.12 

255.39 


158.0 

256.14 

255.89 

255.99 

256.07 

255.63 

255.97 


237.6 

256.86 

256.41 

256.65 

256.81 

256.42 

256.73 


316.8 

257.58 

256.97 

257.36 

257.54 

257.14 

257.49 


396.0 

258.27 

257.60 

258.07 

258.27 

257.95 

258.21 


475.2 

258.92 

258.21 

258.74 

258.94 

258.61 

258.89 


554.4 

259.53 

258.82 

259.37 

259.58 

259.43 

259.53 


792.0 

261.16 

260.45 

261.10 

261.32 

261.07 

261.25 


1029.6 

262.53 

262.18 

262.58 

262.84 

262.72 

262.72 

0.52 

0.0 

255.38 

255.38 

255.38 

255.38 

255.38 

255.38 


158.0 

255.38 

255.39 

255.38 

255.38 

255.30 

255.36 


237.6 

255.39 

255.41 

255.39 

255.38 

255.09 

255.34 


316.8 

255.40 

255.45 

255.41 

255.40 

254.88 

255.32 


396.0 

255.44 

255.50 

255.45 

255.44 

254.98 

255.33 


475.2 

255.51 

255.56 

255.51 

255.50 

255.02 

255.37 


554.4 

255.60 

255.64 

255.60 

255.59 

255.07 

255.46 


792.0 

256.00 

256.02 

255.96 

255.99 

255.41 

255.85 


1029.6 

256.52 

256.34 

256.46 

256.53 

255.99 

256.38 i 
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*■— — “ —■ Series Solution (Ref 31) 

O 16 Element, Full [C] 

/\ 32 Element, Condensed [C] 

Q- - — 16 Element, Condensed [C] 



Figure 4. Accuracy Comparison of Computed Transient Temperature Distributions in Thick 
Axisymmetric Cylinder As Function of [C] Form and Discretization 
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Figure 5. Comparison of Surface Temperature and Integration Step-Size for Different Integration 
Methods and Matrix Forms For Thirty-Two Element Cylinder Discretization 
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Method 2 marched at a consisderably smaller integration step-size (as expected, see Table 4), and con- 
sequently took 50% additional CPU time to achieve the solution end point. Switching to the condensed 
[C] matrix, e = 0.001 , maintained comparable solution accuracy, and achieved the solution in half the 
CPU time of the full [C] solution. This economy results directly from the increased stability of the con- 
densed [C] differential equation system, Eq (2.21), as evidenced by the significantly larger step sizes used 
by QKNINT for the same Method. This speed factor, coupled with the accuracy determination previously 
discussed, indicates that use of condensed [C] is uniformly preferable to full [C] for practical problems. 

Decreasing the accuracy parameter to e = 0.01 for Method 2, and using the condensed [C] form, 
produced the illustrated results which were obtained on the fringe of instability. Note the surface temp- 
erature excursion at t = 1440 sec. (0.4 hr), the immediate reduction in integration step size, and the return 
to an accurate solution. This is a completely automatic operation within QKNINT and requires no user 
insight or control, other than specification of the error parameter. 

■ j ' , ' 

The results of these as well as many other tests indicate that specifying e = 0.01 will in generaly pro- 
duce the most Speedy solutions while maintaining generally adequate accuracy. However, on rare occasions, 
e = 0.01 has allowed errors to propogate beyond return, with the resultant solution diverging to nonsense. 
Setting e = 0.001 has uniformly resulted in smooth stable solution behavior, usually for somewhat increased 
solution CPU times. The specific choice of the option belongs to the user. 

The results shown in Figure 6 additionally illustrate the automatic features of QKNINT, as well 
as indicating the attainment of smooth, well behaved solutions for e = 0.01 . Shown are transient surface 
temperature histories for the thick cylinder problem with different thermal conductivities. The thermal 
conductivity matrix, the counterpart of the stiffness matrix in elasticity, forms the major contribution to 
the Jacobian of the solution differential equation system, Eq (2.21). Hence, the,stability character is 
strongly influenced by changes in thermal conductivity. QKNINT automatically determined a lower- 
bound estimate to the maximum allowable initial integration step size, HMIN, as a function of k as 
tabulated. It similarly obtained a lower-bound estimate of the maximum allowable step-size, HMAX, as _ 
a function of the generated solution. Note that the k = 69.2 w/m°K (40 Btu/hr ft°R) solution found 
the relatively largest HMIN and smallest HMAX; as a consequence, it required the largest CPU execution 
time. This illustrates the important feature of the QKNINT package of completely automatic operation, 
and a built-in capability to tailor the solution constraints to the particular differential equation system 
being solved, as a function of the error control parameter epsilon. 

B. ACCURACY AND SOLUTION SPEED FOR A PRACTICAL ENGINEERING PROBLEM 

The problem selected for illustrating the engineering usefulness of COMOC is the pulsed firing of 
a small, axisymmetric thrust-vector-control rocket motor. Such devices, commonplace in many aerospace, 
missile and spacecraft systems, are subjected to severe cyclic thermal loading from the contained com- 
bustion of propellants. Hence, very large (radial) temperature gradients can be induced, and the structural 
integrity of the nozzle throat is of fundamental importance. Therefore, accuracy of the computed solution 
especially in the throat region, is critically important. i 

The assessment of computational accuracy was made by comparing COMOC solutions to those 
generated by an independent, state-of-the-art, multi-dimensional heat transfer computer program. The 
theoretical foundation of the solution algorithm for this comparison program is the “electric analog” to 
heat conduction, wherein heat capacities are assigned to nodes, and equivalent thermal “resistances” 
between nodes are established based upon the material thermophysical properties and the node distribution. 
An entirely equivalent theoretical basis exists within the MWR finite element theory established in Section 
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1.73 W/m K 
17.30 W/m °K 
69.18 W/m °K 





= 63.05 W/m 2 
= 367° K 


Solution Time (sec) 


Figure 6. Transient Surface Temperature of Thick Axisymmetric Cylinder as Function of Thermal 
Conductivity, Epsilon = 0.0 1 , Method 2, Condensed Matrix 





II, which leads to an apparent “finite difference” formulation. (32) To form this basis, assume that the 
temperature approximation function within a given solution element, T m *, Eq (2.18), can be represented 
by a spatially uniform temperature located at coordinates (£j) within the element, i.e. , 0 m (£j, t). Employing 
a scalar constant for the weighting function, W k , Eq (2.1 6) in conjunction with 0 m (£j,t) yields an 
algorithm form identical in appearance to Eq (2.21). However, to complete the algebra, and since f? m ($j,t) 
is spatially constant, a finite difference approach must be used to evaluate the element thermal conductivity 
matrix, [K1 ] m , Eq (2.20). Typically, a first order difference is used, with the assumption that the co- 
ordinates (£j) are identical with the element geometric centroid (Xj) for interior elements, and that elements 
at a convention boundary are half-thickness and (£j) lies at the surface centroid. 

Figure 7 is an illustration of the subject rocket motor. The overall length of the device is about 
0.15 meter (6 inches); the downstream half was selected as the computational problem domain, span- 
ning from the center of the combustion chamber through the terminus of the nozzle. Also shown in 
Figure 7 is the 104 finite element discretization used for the majority of COMOC computations. This 
discretization was formed from the parent “finite difference” discretization, which employed quadrilateral 
shaped solution domains. (Note the specification of half-width quadrilaterals at the internal and external 
cylindrical convection surfaces.) Each of these domains were bisected to form the triangular element shape 
required for the finite element solution. Hence, the discretizations were approximately uniform for the 
comparison methods. However, as a consequence of the difference between the formalisms, the 52 finite 
difference quadrilaterals yielded 52 node point temperatures for computation, while the 104 finite elements 
yielded 70 node temperatures for integration. 

The boundary condition specifications were vanishing normal temperature gradients at the left and 
right ends of the solution domain (see Figure 7) and an adiabatic outer surface. In actual practice, the inner 
surface of the motor is loaded by radiation from the contained hot gas. While this is capable of solution, 

Eq (2.1 1), the considerable difference in computed surface temperatures between the methods, to be dis- 
cussed, would deleteriously affect the desired comparisons. Hence, for computational purposes, the nozzle 
was loaded convectively by a uniform internal gas temperature of 361 1°K (6040°F) and the convection 
heat transfer coefficient distribution over the surface of interior elements illustrated in Figure 7. The 
maximum energy flux, occurring at the nozzle throat, was about 9.5 x 10 s w/m 2 (3 x 10 6 Btu/hr ft 2 ), 
and this produced throat temperatures at the end of a 1 0-second firing that were comparable with the 
results of laboratory experimentation. The rocket material was selected as a stainless steel, the thermo- 
physical properties of which vary approximately by a factor of two for the temperature range encountered 
during the ten-second firing. 

Illustrated in Figure 8 are the comparisons between the COMOC computed and “finite difference” 
program results for axial temperature distributions on and near the interior and exterior surfaces of the 
rocket motor after a 10-second firing. The node rows for which the COMOC results are plotted are shown 
in Figure 7. The “finite difference” comparison solutions, from theoretical considerations and near 
equilibrium, must be bounded by the respective finite element solution; in addition, they should lie closer 
to the extreme COMOC computed temperatures. The comparison results confirm theoretical predictions 
excellently at both rocket surfaces. The approach to actual surface temperature by the “finite difference” 
solution is poorest where the temperature levels and gradients (hence departure from equilibrium) are largest. 
At the throat, the COMOC results are 1 1 1°K (200° R) higher, and the “finite difference” solution is effectively 
the mean of the actual temperature distribution within the computational cell. The excursion of the “finite 
difference” solution beyond the finite element bound, at the extreme downstream end of the rocket, is 
probably associated with difficulty in establishing the correct approximation to the vanishing gradient 
boundary condition for the skewed quadrilateral zones employed there (Ref. Figure 7). The correct 
bounded behavior occurs at the left end for a similar boundary condition and regular quadrilateral cells. 
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Figure 7. Discretization of Axisymmetric Rocket Motor Into 104 Triangular Cross-Section 

Ring Finite Elements 







Figure 9 illustrates the comparison between the transient temperature distribution at the nozzle 
throat, and at the exterior wall along Node Lineot, Figure 7, as computed by COMOC and the comparison 
computer program. Along the exterior surface, the finite difference solution is continually bounded by 
the finite element solution. The behavior at the nozzle throat is essentially identical except for a modest 
over-shoot in the neighborhood of t = ~ 3 seconds (0.001 hours). This excursion probably results from 
the gross underprediction of surface temperature by the finite difference program during the initial time 
of firing ( over 333°K (600°R) at time equals one second). The most significant departure from equilibrium 
occurs simultaneously with this underprediction. The assumption that the temperature node location (£j) 
of the finite difference analogy occurs at the surface of half-width quadrilateral elements during rapid transients 
is thus observed to be a poor approximation. 

The results of these comparison studies using the two different computer programs for the 
identical problem have verified the theoretical predictions concerning the two different solution form- 
alisms. They show the significant improvement in accuracy that can result from application of finite 
element techniques to problems involving severe transient phenomena. In combination with the accuracy 
and convergence results for the thick cylinder problem, the mathematical theory and computational 
practice for the finite element solution algorithm appears well founded and validated. 

Extensive experimentation has been conducted on the rocket motor problem, using the COMOC 
computer program to assess solution speed and accuracy. Figure 1.0 shows an actual COMOC output 
page showing the temperature (in common units) in the nozzle, after one second of firing, in a format that 
is geometrically similar to the actual rocket motor geometry. The boundary of the motor has been sketched 
in, as well as the approximate location of the 333°K, 61 1°K (600°R, 1 100°R) isotherms. Such user- 
oriented routines greatly enhance the readability and practical usefulness of computer output; these 
routines may be easily prepared by the engineer using COMOC, and can present geometrically similar 
data outputs for any arbitrarily irregular discretization. The user may output the temperature derivative 
vector, |0(t)}’ , Eq (2.21), as well as is illustrated. A very quick glance at this array identifies locations 
where the solution field is most rapidly changing. 

The 104 finite element discretization of the rocket motor is rather coarse. To check accuracy of 
COMOC, a 538 finite element discretization, shown in Figure 1 1 , was established by quartering each 
finite element of the original discretization and adding additional elements in the thick section behind the 
nozzle throat. Shown in Figure 1 2 is the COMOC computed temperature distribution within the motor 
after one-second of firing. The rocket boundary has been sketched in, as well as the 333°K, 61 1°K 
(600°R, 1 100°R) isotherms for comparison to Figure 10. Note also how the output routine has been 
altered to encompass the larger problem output. 

Assessment of the solution accuracy for the 104 element case, made by comparing Figures 10 and 
12, shows that the coarser discretization tended to uniformly underpredict temperature at any given 
location, in agreement with the thick cylinder discretization results. Shown in Figure 13 are the computed 
radial temperature distributions at the nozzle throat and through the combustion chamber wall (Node 
Lines a and 0, Figure 7, respectively). At the nozzle, the 104- element case predictions are 1 1 1°K-167°K 
(200-300°R) lower, although the computed temperature gradients are essentially identical. The 538-element 
case has added significant temperature definition to the region interior of the throat, Node Line a. Very 
little difference in temperature predictions are noted through the combustion chamber wall, Node Line (3. 
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Solution Time - Seconds 


Figure 9. Comparison of Computed Transient Temperature Distribution at Nozzle Throat of Rocket 

Motor During Ten Second Firing 
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Figure 10. Computed Temperature Distribution within Rocket Motor After One Second 
Firing, 104 Finite Element Discretization 
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Figure 1 1. Discretization of Axisymmetric Rocket Motor Into 538 Triangular Cross-Section 

Ring Finite Elements 
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Figure 13. Computed Radial Temperature Distributions Through Rocket Nozzle Throat and 

Combustion Chamber After One Second Firing 





Solution accuracy is affected by the integration method employed as well as discretization. The 
QKNINT package was adapted to allow the Euler integration algorithm (one-step, first-order accurate, 
predictor - no corrector), so favored by practicing engineers, to operate within the automatic error con- 
trol features of QKNINT. Theoretical considerations 00) indicate that both of the QKNINT methods 
enjoy greater accuracy, as well as increased stability interval, as compared to the Euler method. Shown 
in Figure 14 is a comparison of the computed transient temperature solutions at the nozzle throat, for 
the 538 finite element discretization, during the first half second of firing. The smooth curve is faired 
through the QKNINT, Method 1 (QKNINT(l) ) solution. The results of the Euler integration are 
observed to selectively lie above and below the QKNINT values to the maximum of ±5.5°K (±10°R). The 
differences are certainly small, but the QKNINT solution does appear smoother than the Euler predictions. 

Evidence confirming that the QKNINT(l) solution can be more accurate is illustrated in Figure 
15, which presents the computed nozzle surface temperature transient distribution, and the integration 
step size history for the two different integration methods and discretizations. The indicated surface 
temperature is faired through the square symbols, obtained with QKNINT(l) operating in the completely 
automatic mode, for the error parameter, e = 0.01 , and for the integration step size history shown. The 
surface temperature distribution computed by QKNINT(l), starting with an initial step size three orders 
of magnitude larger than the automatic mode start, as illustrated by the circles, is essentially identical 
with the automatic mode results, and was obtained for the step size history illustrated. Therefore, the 
QKNINT(l) automatic mode results were obtained for integration step sizes well within the stability 
bounds of the integration algorithm. Integrating at step sizes nearer, but still within (by epsilon control) 
the bounds, did not degrade the solution. 

The results of the Euler solution to the same problem, shown in Figure 15, indicate surface temper- 
ature predictions uniformly lower than the QKNINT (1) by up to 14 K (25 R). Concurrently, the Euler 
solution predicted somewhat higher solution temperatures along the exterior surface of the rocket motor, 
as illustrated by comparing Figures 16 and 17, which are computer output at the nine-second mark of 
the rocket firing. The step size history for the Euler test is also illustrated in Figure 15. Decreasing the 
error parameter to e = 0.001 did not measurably affect the predicted temperature history of either the 
Euler or QKNINT(l) solutions. Hence, from the proven convergence with discretization, and the results 
of the 538 element test case, the QKNINT(l) solution appears to enjoy greater accuracy than the Euler 
solution, in agreement with theory. 

Shown for illustration purposes on Figure 15 is the transient temperature solution, through one 
second of the rocket motor problem, obtained for the 538 element discretization and QKNINT(l). It 
lies uniformly above the 104-element solution. The plotted step size history confirms the decreased 
stability of the differential equation system, Eq (2.21), resulting from use of a considerably finer solution 
domain discretization. Hence, the improved solution accuracy at the nozzle throat is obtained at a con- 
siderable increase in solution CPU time. 

Shown in Table 6 are comparisons of computer CPU time for the standard 1 0-second firing of the 
104-element rocket motor, as a function of integration algorithm and operation mode of QKNINT. The 
first entry, the standard for comparison, was obtained for completely automatic operation of QKNINT(l). 
The CPU time of 21 .7 seconds was obtained for the case where the user provides no knowledge concerning 
the solution behavior to the computer program. The automatic determination of the minimum initial 
acceptable step size (HMIN) by QKNINT is conservative; inputting a larger HMIN, determined to be 
acceptable by experimentation, decreased solution time to 14.4 second CPU for QKNINT still operating 
under error control (defined as semiautomatic). This solution speed is comparable to the 13.9 sec CPU 
taken by Euler to integrate the same problem in the identical semiautomatic mode, i.e., input HMIN and 
operate under error control. 
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Figure 14. Comparison of Computed Transient Temperature Distribution at the Nozzle Throat for 
QKNINT Method 1 and Euler, 538 Finite Element Discretization 
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Symbol 

Integration 

Method 
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QKNINT (1) 

104 

□ 

QKNINT (1) 

104 

A 

Euler 

104 

V 

QKNINT (1) 

538 



Solution Time (seconds) 


Figure 15. Comparison of Surface Temperature Distribution and Integration Step-Size History for 

QKNINT( 1 ) and Euler, 1 04 and 538 Finite Element 
Discretizations, Epsilon = 0.01 
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QKNINT 

CUMULATIVE NUMBER OF TIMES STEP SIZE HAS DECREASED - 0 AKO INCREASED 

Figure 16. Computed Temperature Distribution within Rocket Motor After Nine Second 
Firing Using Euler Integration Algorithm 
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TABLE 6 

COMPARISON OF SOLUTION SPEED FOR QKNINT (1) AND EULER, 

EPSILON = 0.01 


Symbol 

Figure 

15 

No. of 
Elements 

Integration 

Method 

Operating 

Mode 

CPU Time 
Seconds 

Integration Step-Size 
Extremum, sec. 

TFINAL 

sec. 

H min 

H max 

□ 

104 

QKNINT (1) 

Automatic 

21.7 

0.0596 

1.092 

10.4 

o 

104 

QKNINT (1) 

Semiautomatic 

14.4 

0.360 


10.8 

A 

104 

EULER 

Semiautomatic 

13.9 

0.0144 


10.1 


104 

QKNINT (1) 

Manual 

11.6 

0.360 

1.440 

10.8 


104 

QKNINT (1) 

Manual 

8.9 

0.360 


1.08 


104 

QKNINT (1) 

Manual 

22.7 

0.360 

0.360 

10.4 


104 

QKNINT (1) 

Automatic 

47.6* 

0.0144 

0.461 

10.4 


532 

QKNINT (1) 

Automatic 

162.0 

0.0036 


01.08 


'Epsilon = 0.001 
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The semiautomatic QKNINT (1 ) run established lower bounds on acceptable step size history. Using 
the semiautomatic step-size distribution, and going to completely manual operation (no error control), re- 
duced CPU time to 1 1.6 sec, a 19% reduction over QKNINT (1) semiautomatic and 16% less than Euler. A 
36% reduction over Euler to 8.9 sec CPU was accomplished by manual QKNINT (1 ) operating uniformly at 
HMAXof the semiautomatic mode, except for the first integration step taken at HMIN. Solution accuracy 
was essentially identical to semiautomatic QKNINT (1). As illustration that testing in automatic QKNINT 
(1 ) is not optimal, a run was made at constant step size equal to input HMIN. The CPU time comparison of 
22.7 seconds is just measurably larger than the completely automatic case which attained a maximum step 
size (HMAX) four times larger than HMIN. Decreasing epsilon to 0.001 increased CPU time to 47.6 seconds 
for QKNINT (1) automatic mode, and did not measurably affect the solution. The one-second firing of the 
538 element discretization required 162 seconds CPU for QKNINT (1 ) automatic. Finally, the CPU time for 
the QKNINT (1) fully automatic, 104-finite-element discretization solution was a factor of five smaller than 
the comparison “finite difference” computer program, which employs an Euler integration procedure. 

The solution speed for the completely automatic operation mode of QKNINT was found to depend 
quite strongly on the internal algebraic factor by which current step-size is multiplied during the seeking 
process for HMIN. Shown in Figure 1 8 is the comparison of CPU time to process the standard rocket 
motor test case as a function of the HMIN seeking factor. The optimum value for the factor appears to 
be in the range 3-1/2 to 4. On either side, a too conservative estimate of HMIN is found with a resultant 
increase in CPU time for execution. The automatic mode in QKNINT presently employs a factor of 4. 

It has similarly been determined that requesting output at exact time stations, which requires internal 
interpolation, can increase CPU time by 10-20% dependent upon the problem. 

COMOC operates independent of whether the problem is two-dimensional or axisymmetric. 

Shown in Figure 1 9 is the computed temperature distribution for the rocket motor initial and boundary 
conditions and geometry, except for a two-dimensional problem specification, after 9 seconds of solution 
time. Comparing the results to Figure 17 shows that overall higher temperatures are computed; however, 
the stability character of the two-dimensional equations is essentially identical to the axisymmetric pro- 
blem, with both automatically establishing an identical HMIN and requiring the same CPU time for 
execution. 

C. VERSATILITY ILLUSTRATION FOR COMOC 

The COMOC heat transfer variant enjoys considerable versatility for diverse problem categories in 
addition to the basic boundary condition generality. Inclusive features are element specifiable internal 
heat generation (or removal), general tensor thermal conductivity, variable element thicknesses in two- 
dimensional problems, and multi-material capability, each with specifiable thermophysical properties. 

The transient temperature history of the rocket motor has been illustrated. A disastrous manu- 
facturing defect, amounting to an internal void located behind the nozzle throat, was simulated by 
identifying a new material with zero thermal conductivity in the immediate throat neighborhood elements 
located between Node Rows B-B', Figure 7. The computed solution remained stable in automatic mode 
operation, and the throat surface temperature indicated burn-out after only three seconds of firing. 

One might intuitively estimate that placing an insulator on the exterior rocket surface would raise 
nozzle surface temperatures. The thermophysical properties of elements between Node Rows C-D, Figure 
7, were changed to simulate asbestos and the ten second firing repeated. Interestingly, the nozzle surface 
temperature was not measurably altered, while the exterior surface temperatures (Node Row D) remained 
essentially at the initial 294°K value. Only the temperatures along Node Row C were increased, by a 
maximum of about 75°K. 
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Figure 19. Computed Temperature Distribution in Two-Dimensional Analog of Axisymmetric 

Rocket Motor Problem 



Of significant engineering design value would be evaluation of candidate methods for reducing throat 
temperatures. Of course, using a higher conductivity material with adequate strength properties, if available, 
would accomplish the desired reduction. A more feasible solution would be to install a high thermal 
conductivity insert into the nozzle throat region and keep the remainder of the motor stainless for strength. 
The throat region interior to the small cross-hatched border, Figure 7, was changed to the material pro- 
perties of graphite, and the ten second firing repeated. Shown in Figure 20 are radial temperature dis- 
tributions with and without the insert, at axial stations at the throat (Node Line oc, Figure 7), along the 
graphite-stainless interface upstream of the throat (z = 0.055 m), and within the stainless further upstream 
of the throat (z = 0.051 m),. A definite knee in the temperature is observed at the graphite-stainless inter- 
face on Node Line a, Figure 20. The throat surface temperature is approximately 60°K cooler, at that 
time in the firing, while the material behind the throat is appropriately somewhat hotter. A smaller 
temperature difference is noted, Figure 20, at the upstream graphite-stainless interface. It is noteworthy 
that COMOC computes temperatures at node points lying on a line of material discontinuity with no 
requirement for any user guidance. The final temperature distribution in Figure 20 verifies that the extra 
heat removed by the graphite insert at the throat is conducted throughout the remainder of the motor. 

A promising variation is the hypothesis that the motor be made from material in every regard 
stainless except for its radial thermal conductivity, which approximates that of graphite, i.e., material 
possessing a favorable thermal conductivity tensor. Shown in Figure 21 is the comparison of the axial 
temperature distributions after a ten-second firing for the stainless and pseudo-stainless orthotropic con- 
ductivity rocket motor. A dramatic reduction in interior surface temperature is computed everywhere 
upstream of the throat at the expense of significantly increased temperature levels at the nozzle terminus. 
These results indicate that if such a material were available it might be ideal for rocket motor construction. 

Active cooling, by the passage of propellant through drilled passages near the interior surface of 
the rocket motor, is a proven method for reducing temperature levels in the nozzle throat region. Such 
a problem description is truly three-dimensional; however, a two-dimensional simulation of cooling could 
yield informative design guidance. As an illustration, an extra row of nodes was placed parallel to and 
just exterior, to Node Row B, Figure 7. Additional discretization was added to form a band of triangular 
elements directly behind those on the interior surface, and the presence of active fluid cooling simulated 
by applying negative heat addition therein. Shown in Table 7 is the comparison of computed radial temp- 
erature distributions along Node Line a, Figure 7, resulting from various levels of simulated active cooling 
of the rocket motor during the 1 0 second firing. 

The ability of COMOC to accept a finite element thickness distribution for two-dimensional 
problems renders it applicable to problems where three-dimensional effects exert a prominent influence 
but for which a full three-dimensional solution is not warranted. As an illustration, the 104 element rocket 
motor discretization was considered as the plane cross section of material that had a 2: 1 thickness variation 
over its length (the z axis, Figure 7). Welding was assumed to have occurred between Node Rows A and 
B, Figure 7, within the region 0.046 < z < 0.061 m. The surface temperature distribution occurring 10 
seconds after cool down was initiated is shown in Figure 22.The smaller slope of temperature in the 
thicker region (to the right of the weld centerline) is quite apparent from the dashed mirror image over- 
lay. 
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Figure 20. Radial Temperature Distribution Comparison for Graphite Insert and Stainless Steel 

Rocket Motor 


Temperature - ( K) 
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TABLE 7 


COMPUTED RADIAL TEMPERATURE DISTRIBUTION 
ALONG NODE LINE a, FIGURE 7, FOR SIMULATED ACTIVE COOLING 
OF ROCKET MOTOR AFTER 10 SECOND FIRING 





Temperature °K ( c 

R) 


Node Sequence 
Node Line a 

Simulated Active Cooling Rates — W/m 2 (BTU/hr ft 2 ) 

0 

3.15 x10 s 

(10 6 ) 

3.15 x 10 6 

(10 7 ) 

1 

1810 

(3255) 

1804 

(3245) 

1770 


2 

1511 

(2720) 

1503 

(2706) 

1460 

(2626) 

3 



999 

(1780) 

940 

(1690) 

4 

637 

(1145) 

616 

(1109) 

587 

(1054) 

5 

453 

(798) 

430 

(774) 

416 

(747) 

6 

420 

(755) 


(729) 

392 









Figure 22. Axial Surface Temperature Distribution 1 0 Seconds After Welding 
Operation on Segment of Wedge Shaped Domain 
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VI. CONCLUSIONS AND RECOMMENDATIONS 


The documented computational experiments with the Thermal Analysis Variant of the COMOC 
computer program system have conclusively demonstrated the considerable versatility, engineering use- 
fulness, computational efficiency, and mathematical soundness of the computational embodiment of a 
finite element solution algorithm for solution of transient problems of energy conservation in stationary, 
multi-material continua. Considerable computational speed is indicated using a highly-stable explicit 
numerical integration procedure, for problems of practical engineering importance involving complex 
geometrical shapes and nontrivial boundary condition specifications. The ease with which output can be 
oriented for quick engineering evaluation has been demonstrated. As a system, COMOC appears to 
effectively communicate to the practicing design engineer the sophisticated numerical mathematics of 
modern day digital computing software systems. 

As a natural consequence of the considerable numerical experimentation performed in the de- 
velopment of COMOC, several areas have emerged in which the performance of the computer program 
system could be improved: 

1 . It is apparent that testing for solution accuracy in QKNINT is costly, and that the error- 
dictated discarding of entire integration steps is a major contributing factor. A recently completed 
theoretical study indicates that the accuracy test could be performed at the end of the prediction step 
rather *than after the second correction step. With such a Change in QKNINT, epsilon becomes actual 
percentage error. More importantly, from the economy standpoint, a too large integration interval is 
determined with one-third the execution time of the presently employed method. This would not 
materially affect CPU time for accepted integrations, but would decrease by up to 40% the CPU time 
required to complete an integration at a reduced step-size. This improvement is experimentally verified, 
as detailed in Table 6, for the QKNINT (1 ) manual mode experiments, which decreased by 19% the CPU 
time of the semi-automatic mode solution while marching at the identical step-size distribution. 

2. It is quite apparent that the high stability methods in QKNINT are capable of generating 
accurate solutions using integration step-sizes that are comparable to, and may actually exceed, the de- 
sired output print intervals. If the solution accuracy of an Euler integration is acceptable, the QKNINT 
package could be readily adapted to automatically switch over to Euler when output requirements con- 
strain the integration step size below the stability limit of the highly stable method. A significant savings 
in CPU execution time could then result for problems that are output-bound while retaining the automatic 
error control features within the QKNINT package. 

3. The natural coordinate formulation for finite element matrices has improved accuracy and 
the speed of matrix computation in comparison to previously employed techniques. The extension of 
COMOC to full three-dimensionality using these functional representations is very straightforward and will 
considerably expand its applicability. At the same time, COMOC could be extended to handle problems 
involving a mix of two-and three-dimensional elements such that complex aerospace geometries could be 
readily and economically analyzed. 

4. Active cooling of aerospace structures is of paramount importance in many applications, e.g., 
laser mirrors and hypersonic aerodynamic panels. The finite element formulation for energy conservation 
is readily extended to include a “fluid cooling element” by enforcing both mass and energy conservation, 
for the additional (fluid) material specification within the solution domain. The directional character of 
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the fluid flow would be automatically enforced by the convective additions to the [K] matrix system. 
Hence, an accurate accounting of heating of the coolant and the resultant temperature distribution with- 
in the containing structure would result. The addition of such a “fluid element” capability within 
COMOC, in either two- and/or three-dimensional space would provide a uniquely powerful analysis tool 
for economical engineering design studies of complex aerodynamic and aerospace hardware systems. 
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